%% Check coverage of CIs & Compute empirical coverage probability

function [predict_se, empirical_cov, cov_p_value] = check_CI(hat_vartheta, rep_val_se2, bar_tau, bar_V_tau, nu, cov_stat_sim)
% INPUT: validation study estimates hat_vartheta, reported validation study se^2 rep_val_se2, 
% posterior moments bar_tau, bar_V_tau, EB estimate nu, simulated sample of coverage statistics cov_stat_sim
% OUTPUT: standard error of the predictive distribution predict_se, empirical coverage probability empirical_cov, p value of coverage prob

    N = size(hat_vartheta,1);

    % Compute the upper and lower bounds
    predict_se = sqrt(rep_val_se2 + nu + bar_V_tau);
    lower_bound = bar_tau - 1.282 * predict_se;
    upper_bound = bar_tau + 1.282 * predict_se;

    % Check whether the reported estimates lie in the predicted CIs
    hit = (hat_vartheta >= lower_bound) & (hat_vartheta <= upper_bound);

    % Compute empirical coverage prob
    empirical_cov = sum(hit)/N;
    
    % Compute p-value for coverage statistics
    cov_stat_obs = (empirical_cov - 0.8)^2;
    cov_p_value = mean(cov_stat_sim >= cov_stat_obs);

end